{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0b2a1f13-07b1-4b72-9c81-4da1ba62aec6",
   "metadata": {},
   "source": [
    "# Quantum Volume\n",
    "\n",
    "[Quantum volume](https://en.wikipedia.org/wiki/Quantum_volume) (QV) is a metric for determining the power of a noisy quantum device.  The QV test is performed for a specified qubit number $N$. If the test is passed, then the device can claim a quantum volume of $2^N$.  In practice, the procedure repeats, until the device reaches a qubit number for which the test fails and its greatest passing score is the device's quantum volume. Though imperfect, the test is a reasonable  approximation of the devices usable processing power. This tutorial will demonstrate how CUDA-Q can be used to perform the quantum volume test.\n",
    "\n",
    "\n",
    "The test consists of the following steps (see figure below):\n",
    "1. A special random circuit is constructed (details below)\n",
    "2. A simulation determines the exact probability distribution of every bitstring and the median probability is determined.\n",
    "3. Every bitstring which has an associated probability greater than the median, is considered a heavy bitstring for that particular circuit.\n",
    "4. The circuit is sampled on the noisy device, and the percent of shots resulting in heavy bitstring are\n",
    "5. The process is repeated many times and the resulted averaged. The test is passed if the average is greater than 2/3.\n",
    "\n",
    "\n",
    "\n",
    "<img src=\"./images/QVprocedure.png\" alt=\"Drawing\" style=\"width: 700px;\"/>\n",
    "\n",
    "the circuits used are square, meaning they have the same number of layers as qubits. Each layer consists of a random permutation of qubits, followed by random SU4 operations applied to n/2 pairs of qubits.  (See the first step in the figure above).  For CUDA-Q implementation, the SU4 gates are decomposed using the KAK decomposition (figure from this [paper](https://arxiv.org/pdf/1811.12926)).\n",
    "\n",
    "\n",
    "<img src=\"./images/kakdecomp.png\" alt=\"Drawing\" style=\"width: 500px;\"/>\n",
    "\n",
    "The cell below specifies a circuit size `n` and two CUDA-Q kernels, one performing an SU4 operation and another building the entire QV circuit. This example is constructed for an even number of qubits for simplicity.\n",
    "\n",
    "The QV kernel concludes with application of a bit flip operation on each qubit. This is not part of the QV circuit, but will be used later to introduce noise to the circuit. Otherwise the test would pass every time!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b748ca48-8c65-4158-8c8a-126bbbb39afe",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudaq\n",
    "import numpy as np\n",
    "\n",
    "# Select an even number\n",
    "n = 4\n",
    "su4_per_circuit = int(n / 2 * n)\n",
    "n_params_in_su4 = 21\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def su4_gate(q0: cudaq.qubit, q1: cudaq.qubit, params: list[float]):\n",
    "    u3(params[0], params[1], params[2], q0)\n",
    "    u3(params[3], params[4], params[5], q1)\n",
    "    x.ctrl(q0, q1)\n",
    "    u3(params[6], params[7], params[8], q0)\n",
    "    u3(params[9], params[10], params[11], q1)\n",
    "    x.ctrl(q1, q0)\n",
    "    u3(params[12], params[13], params[14], q0)\n",
    "    x.ctrl(q0, q1)\n",
    "    u3(params[15], params[16], params[17], q0)\n",
    "    u3(params[18], params[19], params[20], q0)\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def qv(n: int, params: list[float], permutations: list[int]):\n",
    "\n",
    "    reg = cudaq.qvector(n)\n",
    "    param_index = 0\n",
    "\n",
    "    for layer in range(n):\n",
    "        for gate in range(n // 2):\n",
    "            su4_gate(reg[permutations[layer * n + gate * 2]],\n",
    "                     reg[permutations[layer * n + gate * 2 + 1]],\n",
    "                     params[param_index:param_index + 21])\n",
    "            param_index += 21\n",
    "\n",
    "    x(reg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3d3fb4f-ca18-4eea-9b65-609fbc31317c",
   "metadata": {},
   "source": [
    "Each circuit must be random. These function randomly choose parameters and permutations for each circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b139613f-f6ca-48de-a4c2-59aa7d705764",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_random_params() -> list[float]:\n",
    "\n",
    "    params = np.random.uniform(0, 2 * np.pi, n_params_in_su4 * su4_per_circuit)\n",
    "\n",
    "    params_list = params.tolist()\n",
    "\n",
    "    return params_list\n",
    "\n",
    "\n",
    "\n",
    "def generate_random_permutations() -> list[int]:\n",
    "\n",
    "    circuit_permutations = []\n",
    "\n",
    "    for i in range(n):\n",
    "        circuit_permutations.extend(\n",
    "            np.random.permutation(n).astype(np.int64).tolist())\n",
    "\n",
    "    return circuit_permutations\n",
    "\n",
    "\n",
    "parameters = generate_random_params()\n",
    "permutations = generate_random_permutations()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d11fd5a-da31-4645-976d-7d944ba4e18e",
   "metadata": {},
   "source": [
    "This function is an auxillary function used later to convert an integer into a \"big endian\" bitstring. This is used to help determine the heavy bitstrings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ded42cc6-3bf7-4b56-90ca-0e69cbe85dad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_bitstring(integer) -> str:\n",
    "\n",
    "    return bin(integer)[2:].zfill(n)[::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d631da6-97e3-473c-b59f-beb2302fa75e",
   "metadata": {},
   "source": [
    "The `percent_heavy_sampled` function takes the random circuit parameters and permutations and the error rate and returns the percent of heavy bitstrings produced by a noisy circuit sample.  \n",
    "\n",
    "The function first sets up the noise model.  It assumes that each $X$ gate applied at the end of the circuit will fail with some probability denoted by the `error_rate`.\n",
    "\n",
    "Next, the noiseless simulation is performed on a GPU simulated with the `nvidia` backend to obtain the state vector.  The `density-matrix-cpu` backend is used to sample the noisy circuit. \n",
    "\n",
    "The rest of function processes these results to determine the heavy bitstring sample probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e852e3c5-29a0-4588-ac56-a102463aa100",
   "metadata": {},
   "outputs": [],
   "source": [
    "def percent_heavy_sampled(circuit_params,\n",
    "                          layer_permutations,\n",
    "                          error_rate,\n",
    "                          print_output=False) -> float:\n",
    "\n",
    "    # Includes option to print results for a single circuit\n",
    "\n",
    "    # Define a bit flip error applied to all qubits\n",
    "    noise = cudaq.NoiseModel()\n",
    "    bf = cudaq.BitFlipChannel(error_rate)\n",
    "    for i in range(n):\n",
    "        noise.add_channel('x', [i], bf)\n",
    "\n",
    "    # Gets noiseless probability distribution\n",
    "    cudaq.set_target(\"nvidia\")\n",
    "    clean_result = np.array(\n",
    "        cudaq.get_state(qv, n, circuit_params, layer_permutations))\n",
    "\n",
    "    # Performs noisy sampling\n",
    "    cudaq.set_target(\"density-matrix-cpu\")\n",
    "    noisy_result = cudaq.sample(qv,\n",
    "                                n,\n",
    "                                circuit_params,\n",
    "                                layer_permutations,\n",
    "                                noise_model=noise,\n",
    "                                shots_count=1000)\n",
    "\n",
    "    # Converts SV amplitudes to probabilities\n",
    "    probs = clean_result * np.conjugate(clean_result)\n",
    "\n",
    "    # Determines the median value\n",
    "    cutoff = np.median(probs).real\n",
    "\n",
    "    if print_output:\n",
    "        print('The Median for this circuit is:')\n",
    "        print(np.median(probs).real)\n",
    "\n",
    "    # Determines if a bitstring is heavy and saves the bitstring in a list if so.\n",
    "    heavy = []\n",
    "    index = 0\n",
    "    circuit_prob = 0\n",
    "\n",
    "    for outcome_prob in probs:\n",
    "        if outcome_prob.real > cutoff:\n",
    "            heavy.append(make_bitstring(index))\n",
    "            circuit_prob += outcome_prob.real\n",
    "        index += 1\n",
    "\n",
    "    if print_output:\n",
    "\n",
    "        print('The heavy bitstrings for this circuit are')\n",
    "        print(heavy)\n",
    "\n",
    "        print('This circuit has an ideal havy sampling probability of:')\n",
    "        print(circuit_prob)\n",
    "\n",
    "    # Determines percent of noisy sample results that are heavy\n",
    "    prob_heavy_in_noisy = 0\n",
    "    for heavy_bitstring in heavy:\n",
    "        prob_heavy_in_noisy += noisy_result.probability(heavy_bitstring)\n",
    "\n",
    "    if print_output:\n",
    "        print('Percent of time noisy sample returned heavy bitstring')\n",
    "        print(prob_heavy_in_noisy)\n",
    "\n",
    "    # Returns this probability\n",
    "    return prob_heavy_in_noisy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8eb4fff-c031-40d1-b521-7aea8bc1014f",
   "metadata": {},
   "source": [
    "You can test a single circuit below to see if it passes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9cd7d79d-f489-426e-97a7-8664e62799a4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Median for this circuit is:\n",
      "0.04363711\n",
      "The heavy bitstrings for this circuit are\n",
      "['0000', '0100', '0010', '1010', '0101', '1101', '0011', '0111']\n",
      "This circuit has an ideal havy sampling probability of:\n",
      "0.8153219893574715\n",
      "Percent of time noisy sample returned heavy bitstring\n",
      "0.488\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.488"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "percent_heavy_sampled(parameters, permutations, 1, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecd2b646-df3e-4b4c-b9be-851f6bdb6116",
   "metadata": {},
   "source": [
    "The true quantum volume is detemined by repeating the process many times and averaging the results. This function repeatedly applies the `percent _heavy_sampled` function for `n_circuit`number of times and prints if the test is passed and returns the average."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f94562c0-7ceb-4609-9e55-ae8791cd8dc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_qv(n_circuits, circuit_size, prob_of_error) -> float:\n",
    "\n",
    "    n = circuit_size\n",
    "    su4_per_circuit = int(n / 2 * n)\n",
    "    number_of_circuits = n_circuits\n",
    "\n",
    "    counter = 0\n",
    "    circuit_results = []\n",
    "\n",
    "    # Loop over n_circuits\n",
    "    while counter < number_of_circuits:\n",
    "        parameters = generate_random_params()\n",
    "        permutations = generate_random_permutations()\n",
    "        circuit_results.append(\n",
    "            percent_heavy_sampled(parameters,\n",
    "                                  permutations,\n",
    "                                  prob_of_error,\n",
    "                                  print_output=False))\n",
    "\n",
    "        counter += 1\n",
    "\n",
    "    # Average the results\n",
    "    score = sum(circuit_results) / len(circuit_results)\n",
    "\n",
    "    print('The score is:')\n",
    "    print(score)\n",
    "\n",
    "    # Determined if QV test is passed\n",
    "    if score > 2 / 3:\n",
    "        print('passed!')\n",
    "        print('Quantum Volume')\n",
    "        print(2**n)\n",
    "\n",
    "    else:\n",
    "        print('failed QV Test')\n",
    "\n",
    "    return score"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2279de84-bb73-48a9-a33b-3aa2efa66573",
   "metadata": {},
   "source": [
    "Try running the QV procedure for 100 four qubit circuits with a 10% chance of error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3cab12da-7cbc-4ff3-80d8-2a81007942ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The score is:\n",
      "0.7280300000000003\n",
      "passed!\n",
      "Quantum Volume\n",
      "16\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.7280300000000003"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 4\n",
    "calc_qv(100, n, .1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afbc7d0a-2507-42e4-9946-3aa93384c645",
   "metadata": {},
   "source": [
    "an interesting benefit of simulation is the ability to explore how noise might affect the QV results.  In this case, the noise model is trivial, but it is still possible to see a relationship between the probability of error in the $X$ gates and the QV outcome. \n",
    "\n",
    "\n",
    "<img src=\"./images/qvplot.png\" alt=\"Drawing\" style=\"width: 900px;\"/>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
